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Abstract 



The presented article contains a 2D mesh generation routine optimized with the Metropolis algorithm. 
The procedure enables to produce meshes with a prescribed size h of elements. These finite element 
meshes can serve as standard discrete patterns for the Finite Element Method (FEM). Appropriate 
meshes together with the FEM approach constitute an effective tool to deal with differential problems. 
Thus, having them both one can solve the 2D Poisson problem, ft can be done for different domains 
being either of a regular (circle, square) or of a non- regular type. The proposed routine is even capable 
to deal with non-convex shapes. 



1 Introduction 



The variety of problems in physics or engineering is formulated by appropriate differential equations with 
some boundary conditions imposed on the desired unknown function or the set of functions. There exists 
a large literature which demonstrates numerical accuracy of the finite element method to deal with such 
issues pp. Clough seems to be the first who introduced the finite elements to standard computational 
procedures [2] . A further historical development and present-day concepts of finite element analysis are 
widely described in references [TJ|3]. In Sec. 2 of the paper and in its Appendixes AjD the mathematical 
concept of the Finite Element Method is presented. 

In presented article the well-known Laplace and Poisson equations will be examined by means of the 
finite element method applied to an appropriate 'mesh'. The class of physical situations in which we 
meet these equations is really broad. Let's recall such problems like heat conduction, seepage through 
porous media, irrotational flow of ideal fluids, distribution of electrical or magnetic potential, torsion of 
prismatic shafts, lubrication of pad bearings and others Therefore, in physics and engineering arises 
a need of some computational methods that allow to solve accurately such a large variety of physical 
situations. 

The considered method completes the above-mentioned task. Particularly, it refers to a standard discrete 
pattern allowing to find an approximate solution to continuum problem. At the beginning, the continuum 
domain is discretized by dividing it into a finite number of elements which properties must be determined 
from an analysis of the physical problem (e. g. as a result of experiments). These studies on particular 
problem allow to construct so-called the stiffness matrix for each element that, for instance, in elasticity 
comprising material properties like stress - strain relationships [23 [5]. Then the corresponding nodal 
loadsWl associated with elements must be found. 

The construction of accurate elements constitutes the subject of a mesh generation recipe proposed by the 
author within the presented article. In many realistic situations, mesh generation is a time-consuming 
and error-prone process because of various levels of geometrical complexity. Over the years, there were 
developed both semi-automatic and fully automatic mesh generators obtained, respectively, by using the 
mapping methods or, on the contrary, algorithms based on the Delaunay triangulation method [BJ, the 
advancing front method [7] and tree methods [8] . It is worth mentioning that the first attempt to create 
fully automatic mesh generator capable to produce valid finite element meshes over arbitrary domains 
has been made by Zienkiewicz and Phillips [9]. 

The advancing front method (AFM) starts from an initial node distribution formed on a basis of the 
domain boundary, and proceeds through a sequential creation of elements within the domain until its 
whole region is completely covered by them. The presented mesh algorithm takes advantage from the 
AFM method as it is demonstrated in Sec. 3. After a node generation along the domain boundary (Sec. 
3.1), in next steps interior of the domain is discretized by adding internal nodes that are generated at 
the same time together with corresponding elements which is similar to Peraire et al. methodology |10) . 
however, positions of these new nodes are chosen differently according to the manner described in Sec. 
3.2. Further steps improve the quality of the mesh by applying the Delaunay criterion to triangular 
elements (Appendix |E| and by a node shifting based on the Metropolis rule (Sec. 4). 



2 The finite element method 
2.1 The mathematical concept of FEM 

The finite elements method (FEM) is based on the idea of division the whole domain fl into a number 
of finite sized elements or subdomains 51* in order to approximate a continuum problem by a behavior 
of an equivalent assembly of discrete finite elements pp . In the presence of a set of elements tt l the total 
integral over the domain fl is represented by the sum of integrals over individual subdomains Q l 

/ Q r(^,..)*- ? £r(.|,..)*^/ r r(.|....)^- ?J r £ (.i,..) ff 

(1) 

where C(u, denotes a differential operator. The continuum problem is posed by appropriate 

differential equations (e. g. Laplace or Poisson equation) and boundary conditions that are imposed on 
the unknown solution <p. The general procedure of FEM is aimed at finding the approximate solution </> 

A Nodes are mainly situated on the boundaries of elements, however, can also be present in their interior. 
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Figure 1: Figure presents the domain fl and its boundary T. The whole domain f2 could be divided into 
subdomains fi 1 with corresponding line segments T 1 being part of the boundary. The idea of division 
into subdomains (elements) constitutes the main concept of the finite element method. 



given by the expansion: 



PN^cfSNj, (2) 



where Nj are shape functions (basis functions or interpolation functions) ^ and all or the most of 
the parameters <f>j remain unknown. After dividing the domain f2, the shape functions are defined locally 
for elements fl l . A typical hnite element is triangular in shape and thus has three main nodes. It is easy 
to demonstrate that triangular subdomains fit better to the boundary T than others e. g. rectangular 
ones (see Fig. [I]). Among the triangular elements family one can find linear, quadratic and cubic elements 
PQ (see also Appendix [A]). A choice of an appropriate type of subdomains depends on a desired order of 
approximation and thus arises directly from the continuum problem. The higher order of element, the 
better approximation. Each triangular element can be described in terms of its area coordinates L\ , L 2 
and L3. There are general rules that govern the transformation from area to cartesian coordinates 

x = L1X1 + L 2 x 2 + £323 

V = L\V\ + L 2 y 2 + L 3 y 3 <=}■ I L 2 I = I y l y 2 y 3 I \ y \ (3) 





1 = L ± +L 2 +L 3 \ L 3 ) V 1 1 1 / \ 1 

where set of pairs {x\, y\), {x 2 , y 2 ), (#3, y 3 ) represents cartesian nodal coordinates. In turn the area 
coordinates are related to shape functions in a manner that depends on the element order. In further 
analysis only the linear triangular elements will be used. For them, the shape functions are simply the area 
coordinates (see Appendix |A[) . Therefore, each pair of shape functions Nl(x, y),Nf(x, y) for k, I = 1, 2, 3 
could be thought as a natural basis of the triangular element. 

2.2 Integral formulas for elements 



We shall consider the linear expression ( 32 ) derived in the Appendix [Bj 

d 2 d 2 



S4> T (Kj> 



) = i ^ I e dx 2 ~ e dy2j ^ x,y ^ dxdy + j S( t ) P dxd y = ° ( 4 ) 



with the boundary condition <p — 7 on T. In such a simply case of integral-differential problems with a 

d 2 d 2 r . n 

differential operator C = —e-^-r; — Ctt^t, the variable <j> in Eq. (4 1 only consists of one scalar function (j> 

ox oy l 

which is the sought solution, while the constant vector f is represented by the last term in the Eq. Q. 

To find the solution for such a problem means to determine the values of 4>(x, y) in the whole domain 

f2. The values of <p on its boundary T are already prescribed to 7. On the other hand, at the very 

beginning (Eq. ([2])) we have postulated that a function could be approximated by an expansion <j) 

given by means of some basis functions N m (x, y), m = 1, . . . , n (for more details see Appendix [A]). Thus 

another possibility to deal with the Poisson problem is just to start from the functional II (Eq. ( |33[ )) and 

dU - . . 

build a set of Euler equations — = — = where m = 1, . . . , n and <f> m approximates value of the solution 



B thc Einstein summation convention 
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4> calculated at the m-th mesh node. 
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(5) 

an 

and after that we calculate the derivative — = — . Moreover, let's simplify our problem by neglecting the 

d(j) m 

last term in the above-presented equation and imposing cf> = 7 on the boundary T instead. In that 
manner, one obtains the expression 



<9n 



/„ f (e f?) E + « (e f E f 4 + ,E «*.} ** - ( a, 



or in a simplified form 



f /' fdNidN m 8NidN m \ , , , , 



It is worth mentioning that some requirements must be imposed on the shape functions N. Namely, if 
n-th order derivatives occur in any term of L then the shape functions have to be such that their n — 1 
derivatives (pay an attention to the Eq. Q) are continuous and finite. Therefore, generally speaking 
C n -i continuity of shape functions must be preserve. 
In turn, after substituting 

f fdNidN" 1 dN l dN m \ , , 
" ~ / e ( ~7*Z~ ~aZ ' aT ~~ HI - ) dxdy (8) 



/ m = / P N m dxdy (9) 



n 



finally one obtains a set of equations 
or in matrix description 



= ^ K F<t> 1 + f m = for m i 1 = 1 • • • » n ( 10 ) 



^3=K<£ + f = 0. (11) 
d<t> 

It is worth noticing that the matrix K is a symmetric one because of the symmetry in exchange of 
subscripts I and m in Eq. (||. 

Now, we are obliged to employ the division of our domain f2 into a set of subdomains Vt 1 . It gives that 

kt - X>r=E /„, ( m! t y) Sf "2'«h™^ a ~ ) -M, (12) 

/ m = J2f m =J2j ni P(x,y)N im (x,y)dxdy. (13) 



Therefore, after the transformation to / subdomains the expression ( 10 ) becomes 



^ #/ m J ^ + f m j = for i = 1, . . . , J and m = 1, . . . , n 



(14) 



<j> = — (K) 1 f (in matrix notation) 



In fact, the summation in Eq. (14) takes into account only these elements fi 1 which contribute to m-th 
node, however, because of the consistency in notation all elements are included in the sum with the 
exception that those AT* functions for which node m does not occur in i-th element are put equal zero. 
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Figure 2: Figure presents the domain f2 and its boundary T after projection to the polygonal domain. 
It has eight boundary nodes and one central node. Comparing both the initial O and the polygonal 51 
domain one can notice that such a simple projection gives rather rough correspondence between them 
a), however, in some cases it could be a sufficient one i. c. when an integrated function changes very 
slowly in some 5-thick neighbourhood of the boundary V b). 



From now, the whole story is to calculate integrals 

Kt= I e( X) y) (^y)dNy,y) + d Nt^y )d Ny,y) ' 

j ni v dx dx d v d v 

= [ dL x \ 1 d J L 2 e(L 1 ,L 2 ,i 3 )|det J l \ (V N n TT T V T N l " 
Jo Jo 



(15) 



f im = [ p(x,y)N m (x,y)dxdy = [ dL x [ ' dL 2 |det J l \ p{L u L 2 , L 3 )N l ™ (L 1; L 2 , L 3 ) 
Jn { Jo Jo 

where N lm = L lm , L 3 = 1 — L\ — L 2 (see Appendix [a| whereas det J 1 - the jacobian of z-th element, 
T matrix together with V operator in new coordinates are evaluated in Appendix [C] An integration 
over the i-th subdomain f2 8 , which is a triangular element with three nodes, enforces the transformation 
from n-dimensional global interpolation to the local interpolations given by means of N ik (x, y) functions 



where ik = 1,2,3. That is why in Eqs (15) new indices ii,i m appear which further are allowed to take 
three possible values 1,2 and 3 for each element i (the local subspace). 

As a next step, the Gauss quadrature is employed to compute above-written integrals numerically as 
it is described in the Appendix |d| And finally, after incorporating boundary conditions to Eq. (14) by 
inserting appropriate boundary values of <f>, the system of equations can be solved. 



3 Mesh generation 
3.1 Initial mesh 

The domain is a set of points in two-dimensional Euclidean plane 3? 2 (see Fig. [TJ. The initial mesh 
should define the shape of the domain f2 or more precisely its boundary dVL. Let's denote the bd(ll) as T. 
It could form a smooth curve (like a circle) or be a polygon. At the beginning it is necessary to assign the 
starting set of nodes belonging to V . Taking into account polygon it is obvious that the initial mesh must 
consist of its vertices, however, in the case of a smooth curve one can choose the initial mesh differently. 
In the article, the author concentrate on the polygonal domains (see Fig. [2]) that can be formed from a 
smooth curve after placing some initial nodes on its boundary T and connecting them by line segments 
(chords). A role of misplaced boundary nodes will be discussed in Sec. 5. 

Let's start with determining the principal rectangular superdomain as a cartesian product [x m i„, x max ] x 
[ymin,ymax] where Vi £ f! : x > x m i n etc. and the following function meshmit (vertices, radius) where 
the variable vertices determines the number of its sides and the second one gives the radius of its cir- 
cumscribed circle. For instance, one can make use of the Octave GNU project (free open source) [TS] 
and creature both the initial nodes p and the initial triangles t arrays in the case of regular polygon of 
N vertices and lying within a circumscribed circle with a given radius, 
function [p, t] = meshi n j t (N, radius) 

1: phi = [0:2*pi/N:2*pi*(N-l)/N]'; 

2: p = [radius*cos(phi), radius*sin(phi)]; 

4: punter = sum (p,l)/size(p,l); 

5: p = [p; p ccntcl ] ; 
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6: for i = 1:(N-1) 

7: t(i,:) = [i, i+1, N+l]; 

8: end 

9: t = [t; 1, N, N+l]; 
end 

where p ccntcr in line 4 denotes the geometrical center of a figure that is an accurate choice for convex 
cases like regular polygons. However, not only convex type problems are available by presented routine. 
One can set explicitly p table of initial points and compute t table on its basis (lines 6-9). In such a case 
the figure center might require to be shifted, for instance, by the formula given in line 13. That center 
displacement is done in respect to the superdomain center, here set as [0, 0] point. The chosen values 
of weight vector depend on the problem. In Fig. [3j meshes for non-convex figures were obtained with 
weight = [0.25,0.75]. 
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If non-convex figure 

pi = p(sum(p.~2 - repmat(p ccntcr , [size(p,l), l])/2, 2) >= 0, :); 

p2 = p(sum(p.~2 - repmat(p contor , [size(p,l), 1]).~2, 2) < 0, :); 

pccntcr = sum ([ we ight(:, l)*pl; weight (:, 2)*p2], l)/size([pl; p2], 1); 
end 



Following further steps of the algorithm presented in next sections, one can obtain meshes for different 
domains Q (see few examples in Fig. [3| . 

Let's introduce a measure that estimates an element area in respect to the prescribed element area S 
designed by the element size h. The measure Sn — S e i em /S gives a normalized area for each element. 
An estimation of the average deviation from assumed value of the element area provides information of 
mesh quality in the case for their fairly uniform distribution. 








Figure 3: Figure shows four domains Q having different shapes. In brackets, finally established set of 
parameters is written: N p - number of mesh points, Ndivisions ~ number of divisions (according to Sec. 
3.2), Sn ~ a normalized average element area are presented; a) regular polygon - square (258, 8, 1.002); 
b) regular polygon with 16 nodes (376, 6, 1.026) which approximates circular shape well; c) non-regular, 
convex figure (315, 8, 1.01); d) non-regular, semi-convex figure (247, 6, 1.071); and two non-regular, 
non-convex figures e) (245, 7, 0.993) and f) (164, 6, 1.0003) both with weight = [0.25, 0.75]. 



3.2 Adding new nodes to the mesh 

In this section, let's start with the procedure that allows us to add new mesh nodes to the existing ones. 
The initial configuration of the nodes were already defined. It must define well the shape of the divided 
area in aspects explained in the description of the Figure [2] These initial nodes are called the constant 
nodes and are kept immobile through the rest of the algorithm steps. Each triangle could be split up into 
two new triangles by adding a new node to its longest bar. To avoid producing triangles much smaller 
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than denned by the element size h only part of them could be broken up i. e. these for which the triangle 
area is one and half times bigger than A. That condition is set in the algorithm by introducing a new 
control parameter C sp ut. The new node is added in the middle of the triangle longest bar. 




Figure 4: Figure presents a division process of non-regular and circular domains together with their 
boundaries. Pictures a) and c) show meshes with new nodes. Some of them are of the illegal type 
(defined in Sec. 3.2). These nodes constitute starting points for next complementary division that 
transforms such not well-defined elements into the correct ones, see pictures b) and d). 

For each triangle 7fe S O for which C^ plit — 1 
Find its longest bar barf ongest 
Calculate a position of a new node p^id 

If the node is the new one 

Add it to the nodes table 

end 

Update triangles table by replacing the old triangle Tk by two new triangles based on 
end 

It is worth mentioning that presented above algorithm is not quite optimal because some of the new 
nodes could produce triangles with one edge divided by a node resulting from splitting up an adjacent 
triangle. Such triangles are not desirable and are denoted as Tuiegai (see Figures^) and c)). Thus the 
previous procedure needs to be improved. Let's add a few extra steps to it: 

For each triangle Tk € ^ perform checking whether it is of Tui ega i type 

If Tk e Tuie 

Split it up into two new properly defined triangles by connecting so-called illegal node 
with the vertex of Tk lying oppositely to it 
Remove the old Tk triangle 
end 
end 

Figures |4j}) and d) show meshes having only desired elements. 
3.3 The boundary of the domain 

The one of the most important issues to definite is the domain boundary. After determining the boundary 
r by the initial constant nodes (lines 1-18 of the presented below algorithm), the next task is to determine 
which new nodes are lying on boundary line segments Y (as it is visible in Fig. These selection is 
done with a help of the following algorithm: 



G 
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Figure 5: Figure presents the square domain divided into a set of new elements Q l with corresponding 
set of line segments T l being its boundary. A way of finding new nodes constitutes the main point of the 
mesh generation process (see Sec. 3.2) while a selection of nodes is perform according to the algorithm 
from Sec. 3.3 a) nodes a,b,c,d have been classified as boundary nodes whereas b) nodes e,f,g have been 
determined as internal nodes. 



// For an initial node table p (nodes from 1 to N) find all pairs of neighbouring vertices: 

1: pairs = zeros([ ], 2); 

2: for i = 1:(N-1) 

3: pairs = [pairs; i i+1]; 

4: end 

5: pairs = [pairs; N 1]; 

// Connect them by a segment line. If x\ — xi ^ then a function y = ax + b exists and one 

// can find pairs a, b for each such a line segment otherwise a vertical line x = a must be found 

6: diff = p(pairs(:,l), :) - p(pairs(:,2), :); 

7: TOL = l.e-5; 

8: for i = l:size(diff, 1) 

9: if diff(i, 1) > TOL || diff(i, 1) < -TOL 

10: coeff(i, 1) = diff(i, 2)./diff(i, 1); 

11: coeff(i, 2) = p(pairs(i, 1), 2) - p(pairs(i, 1), l).*a(i); 

12: else 

13: coeff(i, 1) = p(pairs(i, 1), 1); 

14: coeff(i, 2) = [ ] or a value out of bounds i. e. the principal rectangular superdomain 

15: end 

16: coeff(i, 3) = min(p(pairs(i,:), 2)); 

17: coeff(i, 4) = max(p(pairs(i,:), 2)); 

18: end 
Establish the table of coefficients a, b once. 
For each new node 

Check whether its coordinates [x, y) fulfill any of y — ax + b equations or x — a where y < y<i 
and y > yi 

If yes classify it as boundary node 
else classify it as internal node 
end 

end 



4 Optimization via the Metropolis method 

Let us define the set of mesh triangles il = {7}, j = 1,2,..., M} and a set T l of triangle mesh elements 
to which a node pi belongs. The closest neighbours C(pi) of the mesh point pi are defined as a subset of 
mesh points pj G V 

\/ Pi 3Tcn: Pi G T C{pi) = { Pj : Pj G T for i ? j}. (16) 

Note, that the closest region is not the same what the Voronoi region [T]. Presented definition is needed 
to proceed with the Metropolis algorithm |14) which will be applied in order to adjust triangle's area 
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to the desired size given by the element size h. In turn, a proper triangulation is the essence of the 
finite element method as it is stated in the Sec. 2. Let us divide the whole problem into two different 
tasks. The first one focuses on finding an optimization for mesh elements being the internal elements 
whereas the second one is developed for so-called the edge elements. They are the elements for which 
one triangle's bar belongs to the boundary T of the domain fi. It is assumed that a proper triangulation 
gives a discrete set of triangles Tj which approximates the domain f2 well. 



a) 



b) 




NODES 


N METROPOLIS 


F. = 0.006 
l 


F. = 0.1 
l 


1 


121 


4 


2 


168 


28 


3 


98 


23 


4 


328 


22 


5 


149 


29 




537 


35 


7 


377 


25 



Figure 6: Figure shows an application of the Metropolis algorithm. Picture a) presents initial positions of 
new nodes just after generating them whereas picture b) shows their positions after node shifts according 
to the procedure described in Sec. 4.1 with the following two values of the force strength Fj: 0.006 and 
0.1 applied to each internal node j — I, ... ,7 and temperature set as T = 0.01. The table presents the 
total number of Metropolis steps that was required to obtain the final result shown in b). 



4.1 The internal elements 

Presented method is based on the following algorithm: 

• Define the element size h and consequently the element area A. 

• Initialize the configuration of triangles and then select the internal nodes Vint 
r} i. e. these nodes does not belong to the domain boundary T. 



{pi :p,£PAp^ 



For each node Pi in Vint find its subdomain defined as a set of triangles % to which the node pi 
belongs. 

Perform the Metropolis approach to every internal node pi within its subdomain CI 1 . The Metropo- 
lis algorithm is adopted in order to adjust an area of each triangle in the node's subdomain to 
prescribed value A by shifting the position of the node pi (Fig. [6] demonstrates robustness of the 
Metropolis approach; compare the node distribution in a) and in b)). That adjustment is govern 
by the following rules: 



Find an area of each triangle Ak (where k = 1,2, 
Tji = Pi — Pj for each pj G f2* connected to node pi 



, , K) in fl l together with the vectors 



Calculate the length of each triangle edge \\fji 
element size h i. e. <$||rji|| = — h 
Calculate the new position of the node p^ ew as 



and its deviation from the designed 



i 



(17) 



where Fj are weights corresponding to magnitude of j-th force applying to node pi . Finally, 
they were set to the constant value F. 
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— Find an area of each triangle A% ew in f2* after shifting pi — > p™ ew 

— Apply an energetic measure £ to a sub- mesh f2*. That quantity could be understand in terms 
of a square deviation of a mesh element area from the prescribed element area A. Therefore, 
in the presented paper the S£ is defined as a sum of a discrepancy between each triangle area 
Ak and A after moving node pi and prior it, respectively 

5£ = Y / (,(A n k ew -A) 2 -(Ak-A) 2 ) k = l,2,...,K. (18) 
k 

If the obtained value of an energetic change is lower than zero the change is accepted. Oth- 
erwise, the Metropolis rule is applied i. c. the following condition is checked 

e- 5E / T > r (19) 

where r is a uniformly distributed random number on the unit interval (0, 1) and T denotes 
temperature. 

— The above-presented algorithm is repeated unless an assumed tolerance will be achieved. 

In order to reach a better convergence of the presented method several other improvements could 
be adopted. For instance, the change in the length of the triangle edge could be an additional 
measure of mesh approximation goodness. That condition will ensure a lack of elongated mesh 
elements i. e. elements with very high ratio of its edge lengths (to see such skinny elements look 
at Fig. 



4.2 The boundary elements 

The Metropolis algorithm applied to boundary nodes slightly differs from the above-described case and 
could be summarize in the following steps: 

• Find all the boundary or edge nodes i. e. nodes for which Pk,edge G T. 

• Find triangles in the closest neighbourhood of the considered pk,edge node. Then calculate an area 
of each triangle Ai^dge- 

• Calculate the force acting on each boundary node except the constant nodes and coming from only 
two boundary nodes connected to it. This imposes the following constrain on the motion of the 
fc-th node in order to keep it in the boundary Y 

where 5||rjfc|| is defined as previously. The force is tangential to the boundary F. 

• Similarly, find an area of each triangle A?ld ge after shifting pk edge — > Pledge according to the force 
F k . 

• Adopt the Metropolis energetic condition to the boundary case i.e. 

6£ = J2((^ge-A) 2 -(A l ,edge~A) 2 ) I = 1, 2, . . . , L (21) 
I 

if e - SE/T > r, (accept) 

otherwise, (reject) 

where T denotes temperature and a random number r € J7(0, 1) as previously. 

The main point of this part is to ensure that the boundary nodes are moved just along the boundary F 
(see Appendix [F|. 



9 



5 Results 

Figure [7] presents the square domain (with the edge length equal y/2) initially divided into a set of new 
elements (the upper picture). Then a mesh generation process can follow two different ways. The first 
of them, denoted as (1), is done after switching off the Delaunay procedure and leads to the uniform 
distribution of 512 identical elements with a normalized triangle area SV = 0.902 (equal 0.0039). On 
the other hand, the second way (2) of creating new elements with help of the Delaunay routine gives 
almost uniform mesh with a normalized average area equal 1.015. Thus employing such an optimization 
pattern returns a result closer to desired one whereas the (1) way is much faster and in that particular 
case also does not make use of the Metropolis procedure. That behavior is caused by a symmetry in 
element and in node distribution, therefore no node shifting is needed. That example should clarify why 
other method than merely finding the geometrical center of each element is required to enhance a mesh 
generation routine. 
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Figure 7: Figure presents two different meshes generated on the basis of a square domain (the upper 
picture). The way denoted as (1) is done without the Delaunay and the Metropolis optimization proce- 
dures giving the uniform distribution of 512 identical elements with Sn = 0.902. The second way (2) 
makes use of the Delaunay and the Metropolis routines giving back almost uniform mesh of 459 elements 
with a normalized average area equal 1.015. 



Mesh optimizations 

Figure [8] presents results obtained by enriching the proposed mesh generator by the Metropolis ap- 
proach. Considered meshes were constructed for two non-regular shapes, one of them is also of non- 
convex type Fig. [8k). As it is clearly seen in Fig. [8] a mesh quality was in that way enhanced. However, 
analysis of computed mesh parameters (SV and N e i em ) is not sufficient to explain such mesh improve- 
ment. Thus, to quantify meshes with non-uniformly distributed elements (see upper cases of Fig. [8]), 
the following measure is put forward S var — mean(\S e iem ~ S\ / S). The better fitting to the prescribed 
element area the smaller value of S var . Computed values of S var vary from 0.22 for not optimized meshes 
to 0.15 after employing the Metropolis rule. Note that for both domains a number of very small elements 
definitely decreased (see also Fig. [6]). Moreover, the Metropolis approach offers wide range of feasible 
mesh manipulations that could be achieved by playing with parameters like the shifting force F, the 
condition of ending optimization (assumed tolerance) and temperature T. The shifting force F can differ 
from node to node or can have the same value for all of them. Furthermore the force strength could 
change after each node division due to the decline in element areas. The higher value of the accepted 
force strength F, the faster the mesh generation routine. However, putting too high or too law value of 
F can enormously increase steps of optimization. 

Figure [9] shows above-discussed meshes enhanced by additional switching on another kind of opti- 
mization i. e. the Delaunay criterion (see Appendix Ej). Both routines were applied in the ordered way 
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Figure 8: Figure presents two different meshes generated by the algorithm with the Metropolis opti- 
mization applied just after new node creation. The obtained mesh characteristics are as they follow: a) 



before S N = 1.0589, N p = 150, N e 



' jv 



1.0421, N p = 231, N e , 



= 258^after S N = 1.0114, N, 
416, after S N = 1.057, N p 



p — 160, N e i em 
227, N elem = 406. 



269, and b) before 



i. e.a mesh reconfiguration by the Delaunay method was added before a node shifting via the Metropolis 
routine. This improves final results in both considered cases. 



a) 




b) 




Figure 9: Figure presents two different meshes generated by the algorithm with both the Delaunay and 
the Metropolis optimizations. The obtained mesh characteristics are as they follow: a) SV = 0.98, N p 
= 168, N e i em = 278, and b) S N = 0.99, N p = 244, N elem = 433. 

To examine mesh stability, let's introduce a small perturbation to a mesh configuration obtained by 
the Metropolis procedure. Applying the Metropolis algorithm leads to the global optimum in element 
distribution for a given number of nodes, thus resulted mesh should have the stable configuration. To 
test this, the Delaunay routine was added one more time just after the Metropolis optimization. The 
highest found changes in mesh quality are depicted in Figure |10| The domains are built with meshes 
having a little bit different parameters than earlier. In other tested cases no mesh reconfiguration has 
been detected. The results show that a small exchange in an element configuration in some cases is able 
to slightly modify the mesh. However, the mesh exchange is hardly seen in the figure |10| Thus, on 
the other hand, that example demonstrates the stability of the proposed algorithm and proves that the 
considered mesh generator can be used with confidence. 

FEM solutions 

Having above-generated meshes one can perform some computations on them. Thus, let's solve 
numerically two examples of 2D Poisson problem and then compare them with their exact solutions. The 
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Figure 10: Figure presents two different meshes generated by the main mesh generator (as depicted in 
previous Figures) but with additional Delaunay routine reconfiguring the old meshes after the Metropolis 
optimization. New mesh characteristics are as they follow: a) SV = 1.0003, N p — 164, N e [ em = 272, and 
b) S N = 0.99, N p = 245, N elem = 432. 



numerical procedure is based on the Finite Element Method already described in Sec. 2. Additionally, 
in that way mesh accuracy will be tested. The first considered differential problem is embedded within 
the rectangular domain [—1 1] x [0 1] (it constitutes the mesh) and has the following form: 



d 2 ^ 

dx 2 



d 2 4> 
dy 2 



= 0, 



(22) 



with the boundary conditions (j> = <j>o for x = — 1 and x = 1, and = otherwise. The solution is given 
by the series 

4(/) 1 cosh(n?rx) . 

^O, y) = V T7 — r sm(mry) (23) 

7r * — ' n cosn(n7r) 

n=l,3,5,... 

with N oo. Figure \TT\ depicts a comparison between an approximation of the analytical solution 
(Eq. ( ]22[ ) ) and a FEM result obtained on the mesh. The mesh were tested for two h values: 0.1 and 
0.06. In the case of a rectangular domain the boundary conditions are fulfilled exactly. Therefore no 
boundary perturbation influences the final result. Analysis of the Fig. |11| confirms a good quality of the 
mesh allowing to solve accurately the problem under consideration. 



a) 



y o.5 



b) 



* FEM 

* series 



******* * 




* 


FEM 


* 


series 




•** * * * 

* ** * * * * 
*** * * 
■* * * I * * 



Figure 11: Figure presents two solutions to the Eq. (22) with the boundary condition <fio set as 1. First 



of them has been obtained by the series ([23]) with N — 200 and presents an approximation of the exact 
solution whereas the second one constitutes a FEM approximation. Computations were performed for 
two different meshes. The maximum of \Acj)\ — 0.018 in the a) case for the element size h = 0.06 whereas 
in the case b) max of |A0| = 0.14 for the bigger element size h = 0.1. 
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Let's look on one more differential problem. The following Poisson equation has been solved both 
numerically and analytically within a circular domain of unit radius 



d 2 (f> d 2 <\> 
dx 2 dy 2 



(24) 



The numerical procedure is again based on the FEM approach. The 
2 — y 2 ) . Figure 



with the boundary condition (f> = 
exact solution is given by the expression (j> — 0.25 
versus a numerical one. Their comparison shows that discrepancy between 



1 



12 



presents the analytical result 
them is less than 0.01. Thus, 



both solutions are in very good quantitative agreement, despite the fact that boundary nodes of the used 
mesh (Fig. [3Jd) do not fulfilled the boundary condition precisely (Fig. 12 1). It is caused by imperfection 
in this approximation of designed circular shape (see also Fig. |3j. On the other hand, the second mesh 
(Fig. [12) 3) has higher element size (h = 0.28) than the first one [h — 0.1) and in that way the boundary 
condition is fulfilled exactly on the T. Summarizing, both meshes suffer some weaknesses but they do 
not influence remarkable solutions. 



a) 



b) 



0.3 
0.2 



-0,1 
1 




y 



s 
2 



01 
1 



***** * »"* 

/*** * $ # 



1 -1 



Figure 12: Figure presents comparison between analytical and numerical solutions obtained to Eq. (24 1 
for two meshes with different clement sizes h. The maximum of discrepancy between both solutions was 
also computed and has the following value max \A<f>\ = 0.0095 in the case a) with the element size h = 0.1 
and in the b) case with the element size assumed equal 0.28 max | Acf>\ = 0.0067. 

All figures presented in the paper were prepared using the Matlab package. 



6 Conclusions 

The proposed mesh generator offers a confident way to creature a pretty uniform mesh built with elements 
having desired areas. Mesh optimizations are done by means of the Metropolis algorithm and the 
Delaunay criterion. Finding that computed measures Sn and S var have pretty good values allows to 
classify the presented mesh generator as the good one. The meshes were also examined in the context 
of their stability to some reconfigurations. It was demonstrated that perturbed meshes hardly differ 
from non-perturbed ones. Moreover, the mesh was tested by solving the Poisson problem on it making 
use of the Finite Element Method. The obtained results are in very good quantitative agreement with 
analytical solutions. This additionally underlines the good quality of the proposed mesh generator as 
well as efficiency of the FEM approach to deal with differential problems. 

A The Lagrange polynomials 

The Lagrange polynomials pk (x) are given by the general formula [TJ lllj 

n 

Pfc(a) = II X ~~ l forfc = l,...,n. (25) 

2=1 
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It is clearly seen from the above-given expression that for x — Xk Pk{xk) — 1 and for x — Xj such that 
j k Pk{ x j) — 0. Between nodes values of Pk{x) vary according to the polynomial order i. e. n— 1 which 
is the order of interpolation. Making use of these polynomials one can represent an arbitrary function 
4>{x) as 

<f>(x) = ^p fc (x)0 fc . (26) 

k 

On the other hand, when the interpolated function <j) depends on two spacial coordinates one can define 
basis polynomials in the form 

p m (x, y) = pu(x, y) = pi(x)pj(y), (27) 

where /, J describe raw and column number for the TO-th node in a rectangular lattice (rows align along x 
and columns along y direction, respectively). And consequently, the set {jp\, . . . ,p m , . . . ,p„} is a basis of 
a n-dimensional functional space because each function p m for m = 1, . . . , n equals 1 at the interpolation 
node (x m ,y m ) and at others. It is easy to demonstrate that such functions are orthogonal |llj . Instead 
of spacial coordinates any other coordinates can be considered. In the case of mesh elements the natural 
coordinates are the area coordinates L defined already in the Sec. The mathematical concept of FEM. 
On that basis the shape functions could be constructed as a composition of these three basis polynomials 
i. e. N m (Li, L2, £3) = p^(Li)p h j{L2)p c K {L^) where the values of a, b and c assign the polynomial order 
in each L^-th coordinate and /, J and K denote the m-th node position in a triangular lattice (i. e. in 
the coordinates £1,^2 and L3, respectively). In the pQ could be found a comprehensive description of 
various elements belonging to the triangular family starting from a linear through quadratic to cubic one. 
For simplicity, in these paper only the linear case is looked on. It explicitly means that shape functions 
Nk = Lk(x, y), where k = 1,2, 3, change between two nodes linearly (see Eq. 

B Variational principles 

We shall now look on the left hand of the Eq. ([!]) i. e. the integral expression II 



n = l n C \ X ^fy---) dn - (28) 

We are aimed at determining the appropriate continuous 4> function for which the first variation <5II 
vanish. If 

6U= k [ -^-H[(j) + kt]]) =0 (29) 



=0 



for any Sip then we can say that the expression II is made to be stationary |12j . The function </> is imbed 
in a family of functions (f> + S4> — (f>(x,y) + Krj(x,y) with the parameter k. The variational requirement 



(Eq. ( 29 1) gives vanishing of the first variation for any arbitrary r\. In the presented article, the variational 
problem is limited to the case in which values of desired function <p at the boundary of the region of 
integration i. e. at the boundary curve T are assumed to be prescribed. Generally, the first variation of 
II has the form 

dli dU dU 

m = M 5(f>+ wJ { ^ ) + W v 5 ^ ) + - m 

and vanishes when 

du „ an n ou „ 

=°>^ =0.oT- =°>--' ( 31 ) 



The condition (31) gives the Euler equations. Moreover, if the functional II is quadratic i.e., if all its 
variables and their derivatives are in the maximum power of 2, then the first variation of II has a standard 
linear form 

an = 64> T (k4> + 1) = o, (32) 



which represents Eq. (30), though, in matrix notation. A vector (f> denotes all variational variables i. e. 
<p and its derivatives as it is written in Eq. ( |3T| ). K denotes stiffness matrix (the FEM nomenclature [1]) 
and f is a constant vector (does not depend on <j>). We are interested in finding solutions to the Poisson 
and the Laplace differential equations under some boundary conditions. These classes of differential 
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problems can be represent in such a general linear form (32). Now, we construct a functional II which 
the first variation gives the Poisson-type equation. Firstly, we define the functional II in the form: 

n=/(iefgV + iefgV + ^l^+/( 7 -W (33) 



2 \dx J 2 \dy 



where dY = \J dx 2 + dy 2 . p, 7 and e can be functions of spacial variables x and y. Secondly, we find the 
first variation of II 

m = ( \^S4> x + e^-S(j)y+p5(jAdxdy+ [ (7 - 4>)5cf>dT (34) 



dx dy 
d 

where 84> x = —Sip. And after integration by parts and taking advantage of the Green's theorem [T| one 
ox 

can simplify the above-written equation to the form 

SIL = J 50 I -e0 - e0 + p\ dxdy + j> eS^dT + J (7 - 4>)6</>dr = 0, (35) 

where — — denotes the normal derivative to the boundary T. The expression within the first integral 
on 

constitutes the Poisson equation 

- - e—. j + p = m il (36) 

ox oy z 



whereas the second term in the Eq. ( 35 ) gives the Neumann boundary condition 

= on T (37) 

on 

and the third one represents the Dirichlet boundary condition 

= 7 on T. (38) 

An important note. The above-presented calculation demonstrates a way in which one can incor- 
porate the boundary conditions of Neumann or Dirichlet type into a variational expression IT. However, 
an appropriately formulated boundary-value problem must include only one kind of b.c. (Neumann or 
Dirichlet b.c.) defined on the whole boundary T or is permitted to mix them but in not self-overlapping 
way. Problems solve in Sec. 5 of the article have only the Dirichlet b.c. 



C Transformation in local L-coordinates 



Let us compute the determinant of the jacobian transformation between the global x,y and a local 
L\,L2,L^ coordinate frame. One notices immediately that the problem is degenerate. That is why, we 
introduce a new coordinate z as a linear combination of L\, L2, L3 i. e. z = L\ + L2 + L3. Note that z 



is not an independent coordinate and has a constant value equal 1. 
([3]) we find the jacobian matrix in the form 



After taking into account relations 



J{L\^L2,L^) = 



_ 9(x,y,z) 



d{L 1 ,L 2 ,L 3 ) 



Ox 

dy 
dL x 
dz 
dLl 



dx 
dL2 

dy 
dL 2 

dz 

d~L 2 



dx 
dT 3 

dy 
dL 3 

dz 
dL~ 3 




Furthermore, we have the relation between the jacobian and an element area 

det J = 2A, 

where A denotes the area of a triangle which is based on vertices (xi,yi), (2:2,2/2), (2:3,2/3)- 
we obtain the coordinates transformation rule 



(39) 



(40) 

And finally, 



dxdy — 2AdLidL2, and L 3 = 1 — L\ — L 2 - 



(41) 
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The relation between the gradient operator V in cartesian and in new coordinates is given by: 



' d d~ 




~dL x d 


dx ' dy 




dx dLi 



dL 2 d dL 3 d dL\ d dL 2 d dL 3 d 



1 

2A 



Odd 
dLi ' dL 2 ' dL 3 



dx dL 2 dx dL 3 ' dy dL\ dy dL 2 dy dL 3 

T 



■ 1 ai 










H 




h ) 





d d d 
dLi ' dL 2 ' dL 3 



(42) 



where L k = (a k x + b k y + c k )/(2A) (k = 1,2,3) and a x = y 2 -y 3 ,6i = x 3 -x 2 ,c x = x 2 y 3 ~ x 3 y 2l the rest 
of coefficients is obtained by cyclic permutation of indices 1, 2 and 3. 



D Numerical integration - the Gauss quadrature 

The l.h.s integral I can be approximated by the Q - point Gauss quadrature [U [TBI H3 HH] 

1= dL x dL 2 \detJ\f(L 1: L 2l L z )^Y^f q (L u L 2: L 3 )W q 
Jo Jo x 



(43) 



where W q denotes the weights for q - points of the numerical integration, and can be found in the Table 
5.3 in [TJ. As it was already said, a set of N k (Li, L 2 , L 3 ) shape functions where k = 1, 2, 3 can be used to 
evaluate each / function in the interpolation series which, for instance, in the highest order 10 - nodal 
cubic triangular element has the following form 



f(L lt L 2 ,L 3 ) = J2 N k(Li, L 2 , L 3 )f k + J2 N k (L u L 2 ,L 3 )Af k + N w AAf 1 



(44) 



k=l 



fc=4 



where f k are nodal values of / function, Af k denote departures from linear interpolation for mid-side 
nodes, and A A/ 10 is departure from both previous orders of approximation for the central nod^] [TJ. 
For linear triangular elements only the first term is important which gives an approximation 



f{L 1 ,L 2 ,L 3 ) = Y J L k f k . 



(45) 



fc=i 



Note, that the r.h.s sum (431 does not include the jacobian | det J\ that should be incorporated by the 
weights W q but it is not (in their values given in Table 5.3 from [TJ). Thus let's add the triangle area to 
the above-recalled formula 

Q 

\detJ\/2j2f q {Lx,L 2 ,L 3 )W q (46) 
and in that way we end up with the final expression for the Q - point Gauss quadrature. 



E The Delaunay triangulation algorithm 

The author would like to remind briefly the main points of the Delaunay triangulation method |13j 
together with their numerical implementation using Octave package [15] . Let V = {pi, i = 1,2,... N} be 
set of points in two-dimensional Euclidean plane 5i 2 . They are called forming points of mesh [13] . Let 
us define the triangle T as a set of three mesh points 

T={t j eVJ= 1,2,3}. (47) 

Using the Delaunay criterion one can generate triangulation where no four points from the set of forming 
points V are co-circular: 

Vp 4 e V A Vl i T \\pi - u\\ > p 2 (48) 

where u is the center of the T triangle and p is its radius. 
The proposed algorithm consists of the following steps: 

c it comes from the local Taylor expansion written for finite differences 



l(i 



Figure 13: Figure shows the main idea of the Delaunay criterion, a) Two triangles (with nodes abc and 
acd, respectively) are not Delaunay triangles, b) after exchange of the edge ac to the edge bd two new 
triangles abd and bed replace the old ones. They are both of the Delaunay type. Circles represent the 
Delaunay zones. 



The triangle's bars are given by the following vectors £12,^13,^23 where U 
for each i ^ j and i, j E {1,2,3} . 



tj , ti 



The cross product of each triangle bars defines a plane. The pseudovector A together with its 
projection on the normal to the plane n-dircction A n arc found 



A = t 12 x t 13 

A n = ho A 



(49) 
(50) 



in order to determine the triangle orientation. If the quantity A n > the triangle orientation is 
clockwise, otherwise is counterclockwise. 

• The determinant of the square matrix T>(T) is built on the basis of the set of triangle's nodes given 
by Eq. p7|) 

/ t\ t\ {t\f + {t\f \ 



V{T) = det 



\ 



''2 

MX 
'3 



V\2 



X 

2 



m? + (if) 2 



(51) 



next the following determinant is calculated in order to find out whether a mesh point p t is outside 
or inside the Delaunay zone (see Fig. 13 1 



V(T)i = det 



( tf 
'2 



t y 
tl 



t" 
'3 



{ttf + {t\f l\ 

(tf) 2 + (t y 3 ) 2 
^^2 



3 

2) T (t3) 

\pf Pi (p?) 2 + (pI) 2 1 J 



(52) 



for each pi e V A p L 4 T ■ 



• If for any point pi its T>(T) < the triangle T is not the Delaunay triangle (see Fig. 13 1). Then 
one need to find other triangles in the closest neighbourhood of the triangle T corresponding to 
the number of Pi inside the Delaunay zone and recursively exchange the bars between T and each 
of them (see Fig. 13 d). 

• Finally, checking whether the new two triangles are the Delaunay triangles takes its place. If so, 
new ones are accepted unless the change is canceled. 
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Now, let us present the main points of the algorithm that could be easily implemented using the Octave 
package or Matlab. 

Set of mesh points V = {pi, i = 1,2,..., N} together with their starting configuration forming initial 
triangular mesh fl = {Tj, j = 1,2, ... , M} where M determines the mesh size. 
while 1 
stop = 1 

for each mesh triangle T from j = 1,2, ... ,M 

ensure the clockwise orientation of the triangle (see above) 
for each mesh point pt not belonging to the triangle T 
calculate its V{T~) k 
if it is lower then 
for the triangle T 

find its neighbouring 'j^ lei 9 hbour triangle where the mesh point p^ belongs to 
exchange the common edge between the triangle T and its neighbour in order to form 
two new triangles 

ensure the clockwise orientation of the triangles (see above) 
calculate the newest T>(T new )k for the new triangle 7~ new 
if T>(T„ ew ) k > V{T) k 

accept the change and update the set of mesh triangles 
put stop = 

else reject the change and return to the previous set of mesh triangles 
end 

end 
end 

end 

end 

if stop = 1 break end 
end 

Algorithm ends up with the new triangular mesh £l ne w 
In order to find an orientation of a triangle T one can check the sign of A n (see Eq. ( 50 )). If it is greater 
than the triangle orientation is clockwise unless counterclockwise. In the latter case, to ensure the 
clockwise orientation one can once flip up and down matrix in Eq. (51) then the triangle orientation 
turns into the opposite one. Obviously, this flipping results in the change of the sign of the matrix 
determinant V(T) -> -£>(T). 



F Boundary nodes — remarks 

For each new boundary node p find its boundary neighbours and save them in the tboundary array. 

• Find two nodes among all neighbouring nodes in tboundary table that are closest to the considered 
p node. Then save them in table c i oses t. This works fine for convex figures. 

• If the figure's shape is not of convex type then the algorithm must be more sophisticated. It requires 
to find two nodes that are aligned with the analyzed p node. For example, one can creature an 
array table containing all possible combinations of that node and any other two nodes from the 
tboundary array. Only one combination from the table should be the correct one i. e. it must 
fulfilled conditions describing one of the boundary line segments. Find and save it in table a u gn . 

Having table c i osest or table a i ign construct the following shifty vector 

shift = repmat(-p, [2, 1]) + p(tablei, :); i = closest, align 

shiftdir = shift. /repmat((diag(shift*shift')). "(0.5), [1, 2]); 
It defines the node shift direction and must go along with one of boundary line segments. 
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